//#include"/export/home/garykche/development/mpi_lasso/src/hmm_impute_dimensions.h"
//#pragma OPENCL EXTENSION cl_amd_fp64 : enable
#include"/home/garykche/gary/code/src/hmm_impute_dimensions.h"
//pragma OPENCL EXTENSION cl_khr_fp64 : enable

#define EMISSION_VEC (local_call_probs[0] * local_penetrance_matrix[s] + local_call_probs[1] * local_penetrance_matrix[3+s] + local_call_probs[2] * local_penetrance_matrix[6+s])

// THIS SECTION IS FOR THE COMPUTE SIGNPOSTS

__kernel void init_forwardbackward_float(
__const int total_template_haplotypes,
__const int markers,
__const int states,
__constant int * current_subject,
__constant int * current_marker_index,
__constant int * offset,
__global float * call_probs,
__global float * penetrance_matrix,
__global float * emission_float_vector,
__global float * current_prob_vector,
__global float * float_block,
__global packedshort_t * short_template_haplotypes,
__global float * haplotype_marginals,
__local float * local_call_probs, // this is a three element array
__local float * local_penetrance_matrix,
__local float * local_haplotype_marginals,
__local packedshort_t * local_short_hap_vec
){
  int hap1 = get_group_id(0);
  if (hap1>=total_template_haplotypes){
    return;
  }
  int threadindex = get_local_id(0);
  //int marker = current_marker_index[0];
  int marker = current_marker_index[0]+1;
  if (threadindex<3){
    local_call_probs[threadindex] = 
    call_probs[*current_subject * 3*markers + 3*marker+threadindex];
  }
  if (threadindex<9){
    local_penetrance_matrix[threadindex] = 
    penetrance_matrix[9*marker+threadindex];
  }
  local_haplotype_marginals[threadindex] = 0;
  int s1 = short_template_haplotypes[(marker*total_template_haplotypes+hap1)/2].compressed[hap1%2];
  barrier(CLK_LOCAL_MEM_FENCE);
  for(int chunk=0;chunk<=total_template_haplotypes/BLOCK_WIDTH;++chunk){
    int hapoffset = chunk*BLOCK_WIDTH;
    int hap2 = hapoffset+threadindex;
    if (hap2<total_template_haplotypes){
      int state = hap1*total_template_haplotypes+hap2;
      int chunkwidth = total_template_haplotypes-hapoffset<BLOCK_WIDTH?total_template_haplotypes-hapoffset:BLOCK_WIDTH;
      if (threadindex<chunkwidth/2){
        local_short_hap_vec[threadindex] =  short_template_haplotypes[(marker*total_template_haplotypes+hapoffset)/2+threadindex];
      }
      barrier(CLK_LOCAL_MEM_FENCE);
      int s = s1+ local_short_hap_vec[threadindex/2].compressed[threadindex%2];
      float conditioned;
      if (offset[0]>-1){
        conditioned =  EMISSION_VEC * float_block[(marker-offset[0])*states+state];
        float_block[(marker-offset[0]-1)*states+state]=conditioned;
      }else{
        conditioned = EMISSION_VEC * current_prob_vector[state];
        current_prob_vector[state] = conditioned;
      }
      // FOR DEBUGGING BELOW
      //emission_float_vector[hap1*total_template_haplotypes+hap2] =  EMISSION_VEC;

      local_haplotype_marginals[threadindex]+=conditioned;
    }
    barrier(CLK_LOCAL_MEM_FENCE);
  }
  for(int s=BLOCK_WIDTH/2; s>0; s>>=1) {
    if (threadindex<s) {
      local_haplotype_marginals[threadindex] += local_haplotype_marginals[threadindex+s];
    }
    barrier(CLK_LOCAL_MEM_FENCE);
  } 
  if (threadindex==0) {
    //haplotype_marginals[hap1] = hap1;
    haplotype_marginals[hap1] = local_haplotype_marginals[threadindex];
    //haplotype_marginals[((marker-1)%2)*total_template_haplotypes+hap1] = local_haplotype_marginals[threadindex];
  }
  return;
}

__kernel void init_forwardbackward_short(
__const int total_template_haplotypes,
__const int total_template_haplotype_pairs,
__const int markers,
__const int compressed_states,
__constant int * current_subject,
__constant int * current_marker_index,
__constant int * offset,
__global float * call_probs,
__global float * penetrance_matrix,
__global float * emission_float_vector,
__global packedshort_t * current_short_vector,
__global packedshort_t * short_block,
__global packedshort_t * short_template_haplotypes,
__global float * haplotype_marginals,
__local float * local_call_probs, // this is a three element array
__local float * local_penetrance_matrix,
__local float * local_haplotype_marginals,
__local packedshort_t * local_short_prob_vec,
__local packedshort_t * local_short_hap_vec
//__local float * local_prob_vec_float,
){
  int hap1 = get_group_id(0);
  if (hap1>=total_template_haplotypes){
    return;
  }
  int threadindex = get_local_id(0);
  int marker = current_marker_index[0];
  if (threadindex<3){
    local_call_probs[threadindex] = 
    call_probs[*current_subject * 3*markers + 3*marker+threadindex];
  }
  if (threadindex<9){
    local_penetrance_matrix[threadindex] = 
    penetrance_matrix[9*marker+threadindex];
  }
  local_haplotype_marginals[threadindex] = 0;
  int s1 = short_template_haplotypes[marker*total_template_haplotype_pairs+(hap1/2)].compressed[hap1%2];
  barrier(CLK_LOCAL_MEM_FENCE);
  for(int chunk=0;chunk<=total_template_haplotypes/BLOCK_WIDTH;++chunk){
    int hapoffset = chunk*BLOCK_WIDTH;
    int hap2 = hapoffset+threadindex;
    if (hap2<total_template_haplotypes){
      int chunkwidth = total_template_haplotypes-hapoffset<BLOCK_WIDTH?total_template_haplotypes-hapoffset:BLOCK_WIDTH;
      if (threadindex<chunkwidth/2){
        local_short_hap_vec[threadindex] =  short_template_haplotypes[marker*total_template_haplotype_pairs+hapoffset/2+threadindex];
        if (offset[0]>-1){
          local_short_prob_vec[threadindex] =  short_block[(marker-offset[0])*compressed_states+hap1*total_template_haplotype_pairs+hapoffset/2+threadindex];
        }else{
          local_short_prob_vec[threadindex] = current_short_vector[hap1*total_template_haplotype_pairs+hapoffset/2+threadindex];
        }
      }
      barrier(CLK_LOCAL_MEM_FENCE);
      int s = s1+ local_short_hap_vec[threadindex/2].compressed[threadindex%2];
      // FOR DEBUGGING BELOW
      //emission_float_vector[hap1*total_template_haplotypes+hap2] =  EMISSION_VEC;
      //return;
      float conditioned = EMISSION_VEC * short2float(local_short_prob_vec[threadindex/2].compressed[threadindex%2]);
      local_haplotype_marginals[threadindex]+=conditioned;
      local_short_prob_vec[threadindex/2].compressed[threadindex%2] = float2short(conditioned);
      barrier(CLK_LOCAL_MEM_FENCE);
      if (threadindex<chunkwidth/2){
        if (offset[0]>-1){
          short_block[(marker-1-offset[0])*compressed_states+hap1*total_template_haplotype_pairs+hapoffset/2+threadindex]=local_short_prob_vec[threadindex];
        }else{
          current_short_vector[hap1*total_template_haplotype_pairs+hapoffset/2+threadindex] = local_short_prob_vec[threadindex];
        }
      }
    }
    barrier(CLK_LOCAL_MEM_FENCE);
  }
  for(int s=BLOCK_WIDTH/2; s>0; s>>=1) {
    if (threadindex<s) {
      local_haplotype_marginals[threadindex] += local_haplotype_marginals[threadindex+s];
    }
    barrier(CLK_LOCAL_MEM_FENCE);
  } 
  if (threadindex==0) {
    haplotype_marginals[((marker-1)%2)*total_template_haplotypes+hap1] = local_haplotype_marginals[threadindex];
  }
  return;
}


__kernel void reduce_fb_marginals(
__const int total_template_haplotypes,
__constant int * current_marker_index,
__global float * haplotype_marginals,
__global float * full_marginal,
__local float * local_haplotype_marginals
){
  int threadindex = get_local_id(0);
  local_haplotype_marginals[threadindex] = 0;
  barrier(CLK_LOCAL_MEM_FENCE);
  for (int chunk = 0;chunk<=total_template_haplotypes/BLOCK_WIDTH;++chunk){
    int haploindex = chunk*BLOCK_WIDTH+threadindex;
    if (haploindex<total_template_haplotypes){
      local_haplotype_marginals[threadindex] += haplotype_marginals[haploindex];
      //local_haplotype_marginals[threadindex] += haplotype_marginals[(current_marker_index[0]%2)*total_template_haplotypes+haploindex];
    }
    barrier(CLK_LOCAL_MEM_FENCE);
  }
  for(int s=BLOCK_WIDTH/2; s>0; s>>=1) {
    if (threadindex<s) {
      local_haplotype_marginals[threadindex]+=local_haplotype_marginals[threadindex+s];
    }
    barrier(CLK_LOCAL_MEM_FENCE);
  }
  if (threadindex==0){
    full_marginal[threadindex] = local_haplotype_marginals[threadindex]; 
  }
  return;
}

__kernel void forward_backward_emission_float(
__const int total_template_haplotypes,
__const int markers,
__const int states,
__constant int * current_subject,
__constant int * current_marker_index,
__constant int * offset,
__constant int * writesignpost,
__global float * site_thetas,
__global int * underflows,
__global float * call_probs,
__global float * penetrance_matrix,
__global packedshort_t * short_template_haplotypes,
__global float * emission_float_vector,
__global float * haplotype_marginals,
__constant float * full_marginal,
__global float * current_float_vector,
__global float * float_block,
__global float * unconditioned_float_vector,
__local float * local_call_probs, // this is a three element array
__local float * local_penetrance_matrix,
__local float * local_haplotype_marginals,
__local float * local_prob_vec,
__local packedshort_t * local_short_hap_vec
){

  bool odd = current_marker_index[0]%2;
  int hap1 = get_group_id(0);
  if (hap1>=total_template_haplotypes){
    return;
  }
  int threadindex = get_local_id(0);
  if (threadindex<3){
    local_call_probs[threadindex] = 
    call_probs[current_subject[0] * 3*markers + 
    3*current_marker_index[0]+threadindex];
  }
  if (threadindex<9){
    local_penetrance_matrix[threadindex] = 
    penetrance_matrix[9*current_marker_index[0]+threadindex];
  }
  local_haplotype_marginals[threadindex] = 0;
  int s = short_template_haplotypes[(current_marker_index[0]*total_template_haplotypes+hap1)/2].compressed[hap1%2];

  float hap1marginal = haplotype_marginals[odd*total_template_haplotypes+hap1];
  float site_theta = site_thetas[*current_marker_index];
  float change0 = (1-site_theta)*(1-site_theta);
  float change1 = (1-site_theta)*(site_theta)/total_template_haplotypes;
  float change2 = *full_marginal * (site_theta)*(site_theta)/(total_template_haplotypes*total_template_haplotypes);
  if (*full_marginal<SMALL){
    change0*=LARGE;
    change1*=LARGE;
    change2*=LARGE;
    if (get_global_id(0)==0) *underflows += 1;
  }
  barrier(CLK_LOCAL_MEM_FENCE);
  for(int chunk=0;chunk<=total_template_haplotypes/BLOCK_WIDTH;++chunk){
    int hapoffset = chunk*BLOCK_WIDTH;
    int hap2 = hapoffset+threadindex;
    if (hap2<total_template_haplotypes){
      int state = hap1*total_template_haplotypes+hap2;
      int chunkwidth = total_template_haplotypes-hapoffset<BLOCK_WIDTH?total_template_haplotypes-hapoffset:BLOCK_WIDTH;
      if (threadindex<chunkwidth/2){
        local_short_hap_vec[threadindex] =  short_template_haplotypes[(current_marker_index[0]*total_template_haplotypes+hapoffset)/2+threadindex];
      }
      barrier(CLK_LOCAL_MEM_FENCE);
      s += local_short_hap_vec[threadindex/2].compressed[threadindex%2];

      float hap2marginal = haplotype_marginals[odd*total_template_haplotypes+hap2];
      if(offset[0]>-1){
        local_prob_vec[threadindex] = hap1==hap2 ? (float_block[(current_marker_index[0]-offset[0])*states+state]  * change0 + change1 * hap1marginal + change2) : (float_block[(current_marker_index[0]-offset[0])*states+state]  * change0 + change1 * (hap1marginal + hap2marginal) + change2 * 2);
      }else{
        local_prob_vec[threadindex] = hap1==hap2 ? (current_float_vector[state] * change0 + change1 * hap1marginal + change2) : (current_float_vector[state] * change0 + change1 * (hap1marginal + hap2marginal) + change2 * 2);
      }
      float conditioned  =  EMISSION_VEC*local_prob_vec[threadindex];
      if (offset[0]>-1){
        // store the final value after including all transition probs
        float_block[(current_marker_index[0]-offset[0])*states+state]=local_prob_vec[threadindex];
        // store the next vector that will include transition probs
        if (current_marker_index[0]-offset[0]>0){
          float_block[(current_marker_index[0]-offset[0]-1)*states+state]=conditioned;
        }
      }else{
        // store the next vector that will include transition probs
        current_float_vector[state] = conditioned;
        if (writesignpost){
          unconditioned_float_vector[state]=local_prob_vec[threadindex];
        }
      }
     // return;
      local_haplotype_marginals[threadindex]+= conditioned;
    }
    barrier(CLK_LOCAL_MEM_FENCE);
  }
  for(int s=BLOCK_WIDTH/2; s>0; s>>=1) {
    if (threadindex<s) {
      local_haplotype_marginals[threadindex] += local_haplotype_marginals[threadindex+s];
    }
    barrier(CLK_LOCAL_MEM_FENCE);
  } 
  if (threadindex==0) {
// modify this to use odd/even
    odd = (current_marker_index[0]-1)%2;
    haplotype_marginals[(odd)*total_template_haplotypes+hap1] = local_haplotype_marginals[threadindex];
  }
  return;
}

__kernel void forward_backward_float(
__const int total_template_haplotypes,
__const int states,
__constant int * current_marker_index,
__constant int * offset,
__global float * site_thetas,
__global int * underflows,
__global float * emission_float_vector,
__global float * haplotype_marginals,
__constant float * full_marginal,
__global float * current_float_vector,
__global float * float_block
){

  int hap1 = get_group_id(0);
  if (hap1>=total_template_haplotypes){
    return;
  }
  int threadindex = get_local_id(0);
  float hap1marginal = haplotype_marginals[hap1];
  float site_theta = site_thetas[*current_marker_index];
  float change0 = (1-site_theta)*(1-site_theta);
  float change1 = (1-site_theta)*(site_theta)/total_template_haplotypes;
  float change2 = *full_marginal * (site_theta)*(site_theta)/(total_template_haplotypes*total_template_haplotypes);
  if (*full_marginal<SMALL){
    change0*=LARGE;
    change1*=LARGE;
    change2*=LARGE;
    if (get_global_id(0)==0) *underflows += 1;
  }
  for(int chunk=0;chunk<=total_template_haplotypes/BLOCK_WIDTH;++chunk){
    int hapoffset = chunk*BLOCK_WIDTH;
    int hap2 = hapoffset+threadindex;
    if (hap2<total_template_haplotypes){
      int state = hap1*total_template_haplotypes+hap2;
      float hap2marginal = haplotype_marginals[hap2];
      if(offset[0]>-1){
        float_block[(current_marker_index[0]-offset[0])*states+state] = hap1==hap2 ? (float_block[(current_marker_index[0]-offset[0])*states+state]  * change0 + change1 * hap1marginal + change2) : (float_block[(current_marker_index[0]-offset[0])*states+state]  * change0 + change1 * (hap1marginal + hap2marginal) + change2 * 2);
      }else{
        current_float_vector[state] = hap1==hap2 ? (current_float_vector[state] * change0 + change1 * hap1marginal + change2) : (current_float_vector[state] * change0 + change1 * (hap1marginal + hap2marginal) + change2 * 2);
      }
    }
  }
  return;
}

__kernel void forward_backward_short(
__const int total_template_haplotypes,
__const int total_template_haplotype_pairs,
__const int compressed_states,
__constant int * current_marker_index,
__constant int * offset,
__global float * site_thetas,
__global int * underflows,
__global float * emission_float_vector,
__global float * haplotype_marginals,
__constant float * full_marginal,
__global packedshort_t * current_short_vector,
__global packedshort_t * short_block,
__local packedshort_t * local_short_prob_vec
){

  int hap1 = get_group_id(0);
  if (hap1>=total_template_haplotypes){
    return;
  }
  int threadindex = get_local_id(0);
  float hap1marginal = haplotype_marginals[hap1];
  float site_theta = site_thetas[*current_marker_index];
  float change0 = (1-site_theta)*(1-site_theta);
  float change1 = (1-site_theta)*(site_theta)/total_template_haplotypes;
  float change2 = *full_marginal * (site_theta)*(site_theta)/(total_template_haplotypes*total_template_haplotypes);
  if (*full_marginal<SMALL){
    change0*=LARGE;
    change1*=LARGE;
    change2*=LARGE;
    if (get_global_id(0)==0) *underflows += 1;
  }
  for(int chunk=0;chunk<=total_template_haplotypes/BLOCK_WIDTH;++chunk){
    int hapoffset = chunk*BLOCK_WIDTH;
    int hap2 = hapoffset+threadindex;
    if (hap2<total_template_haplotypes){
      float hap2marginal = haplotype_marginals[hap2];
      int chunkwidth = total_template_haplotypes-hapoffset<BLOCK_WIDTH?total_template_haplotypes-hapoffset:BLOCK_WIDTH;
      if (threadindex<chunkwidth/2){
        if (offset[0]>-1){
          local_short_prob_vec[threadindex] =  short_block[(current_marker_index[0]-offset[0])*compressed_states+hap1*total_template_haplotype_pairs+hapoffset/2+threadindex];
        }else{
          local_short_prob_vec[threadindex] = current_short_vector[hap1*total_template_haplotype_pairs+hapoffset/2+threadindex];
        }
      }
      barrier(CLK_LOCAL_MEM_FENCE);
      float newprob = hap1==hap2?(short2float(local_short_prob_vec[threadindex/2].compressed[threadindex%2]) * change0 + change1 * hap1marginal + change2):(short2float(local_short_prob_vec[threadindex/2].compressed[threadindex%2]) * change0 + change1 * (hap1marginal + hap2marginal) + change2 * 2);
      local_short_prob_vec[threadindex/2].compressed[threadindex%2] = float2short(newprob);
      // FOR DEBUGGING BELOW
      //emission_float_vector[hap1*total_template_haplotypes+hap2] =  newprob;
      //return;
      barrier(CLK_LOCAL_MEM_FENCE);
      if (threadindex<chunkwidth/2){
        if (offset[0]>-1){
          short_block[(current_marker_index[0]-offset[0])*compressed_states+hap1*total_template_haplotype_pairs+hapoffset/2+threadindex]=local_short_prob_vec[threadindex];
        }else{
          current_short_vector[hap1*total_template_haplotype_pairs+hapoffset/2+threadindex] = local_short_prob_vec[threadindex];
        }
      }
    }
    barrier(CLK_LOCAL_MEM_FENCE);
  }
  return;
}

__kernel void loglikelihood_float(
__const int states,
__constant int * underflows,
__global float * loglikelihood,
__global float * prior_vector,
__global float * block_matrix,
__local float * local_normalizer
){
  int threadindex = get_local_id(0);
  local_normalizer[threadindex] = 0;
  for(int chunk=0;chunk<=states/BLOCK_WIDTH;++chunk){
    int state_offset = chunk*BLOCK_WIDTH;
    int state = state_offset + threadindex;
    if(state<states){
      local_normalizer[threadindex]+=
      block_matrix[state] * prior_vector[state];
    }
    barrier(CLK_LOCAL_MEM_FENCE);
  }
  for(int s=BLOCK_WIDTH/2; s>0; s>>=1) {
    if (threadindex<s) {
      local_normalizer[threadindex]+=local_normalizer[threadindex+s];
    }
    barrier(CLK_LOCAL_MEM_FENCE);
  }
  if (threadindex==0){
    loglikelihood[threadindex] = log(local_normalizer[threadindex])-log(LARGE)*underflows[0];
  }
  return;
}

__kernel void loglikelihood_short(
__const int states,
__constant int * underflows,
__global float * loglikelihood,
__global packedshort_t * short_prior_vector,
__global packedshort_t * short_block_matrix,
__local packedshort_t * local_short1,
__local packedshort_t * local_short2,
__local float * local_normalizer
){
  int threadindex = get_local_id(0);
  local_normalizer[threadindex] = 0;
  for(int chunk=0;chunk<=states/BLOCK_WIDTH;++chunk){
    int state_offset = chunk*BLOCK_WIDTH;
    int state = state_offset + threadindex;
    if(state<states){
      int chunkwidth = states-state_offset<BLOCK_WIDTH?states-state_offset:BLOCK_WIDTH;
      if (threadindex<chunkwidth/2){
        local_short1[threadindex] =  short_prior_vector[state_offset/2+threadindex];
        local_short2[threadindex] =  short_block_matrix[state_offset/2+threadindex];
      }
      barrier(CLK_LOCAL_MEM_FENCE);
      local_normalizer[threadindex]+=
      local_short1[threadindex/2].compressed[threadindex%2] *  
      local_short2[threadindex/2].compressed[threadindex%2];
    }
    barrier(CLK_LOCAL_MEM_FENCE);
  }
  for(int s=BLOCK_WIDTH/2; s>0; s>>=1) {
    if (threadindex<s) {
      local_normalizer[threadindex]+=local_normalizer[threadindex+s];
    }
    barrier(CLK_LOCAL_MEM_FENCE);
  }
  if (threadindex==0){
    loglikelihood[threadindex] = log(local_normalizer[threadindex])-log(LARGE)*underflows[0];
  }
  return;
}


__kernel void collect_xover_marginals_block_float(
__const int states,
__const int total_template_haplotypes,
__constant int * block_marker,
__constant int * hap_pair,
__global float * float_block,
__global float * backward_xover_marginals_big,
__local float * local_prob_vec,
__local float * local_marginals0,
__local float * local_marginals1,
__local float * local_marginals2,
__local float * local_marginals3
){
  int hap1 = get_group_id(0); 
  if (hap1>=total_template_haplotypes) return;
  int state_offset = hap1*total_template_haplotypes;
  int threadindex = get_local_id(0);
  local_marginals0[threadindex] = 0;
  local_marginals1[threadindex] = 0;
  local_marginals2[threadindex] = 0;
  local_marginals3[threadindex] = 0;
  barrier(CLK_LOCAL_MEM_FENCE);
  for(int chunk=0;chunk<=total_template_haplotypes/BLOCK_WIDTH;++chunk){
    int hapoffset = chunk*BLOCK_WIDTH;
    int hap2 = hapoffset+threadindex;
    if (hap2<total_template_haplotypes){
      int state = state_offset + hap2;
      local_prob_vec[threadindex] = float_block[block_marker[0]*states+state];
      barrier(CLK_LOCAL_MEM_FENCE);
      local_marginals0[threadindex] += (hap_pair[HAP_FIRST]==hap1 &&
      hap_pair[HAP_SECOND]==hap2)*local_prob_vec[threadindex];
      local_marginals1[threadindex] += (hap_pair[HAP_FIRST]==hap1) *local_prob_vec[threadindex];
      local_marginals2[threadindex] += (hap_pair[HAP_SECOND]==hap2) *local_prob_vec[threadindex];
      local_marginals3[threadindex] += (hap1<=hap2)* local_prob_vec[threadindex];
    }
    barrier(CLK_LOCAL_MEM_FENCE);
  }
  for(int s=BLOCK_WIDTH/2; s>0; s>>=1) {
    if (threadindex<s) {
      local_marginals0[threadindex]+=local_marginals0[threadindex+s];
      local_marginals1[threadindex]+=local_marginals1[threadindex+s];
      local_marginals2[threadindex]+=local_marginals2[threadindex+s];
      local_marginals3[threadindex]+=local_marginals3[threadindex+s];
    }
    barrier(CLK_LOCAL_MEM_FENCE);
  }
  if(threadindex==0){
    backward_xover_marginals_big[0*total_template_haplotypes+hap1]
    = local_marginals0[threadindex];
    backward_xover_marginals_big[1*total_template_haplotypes+hap1]
    = local_marginals1[threadindex];
    backward_xover_marginals_big[2*total_template_haplotypes+hap1]
    = local_marginals2[threadindex];
    backward_xover_marginals_big[3*total_template_haplotypes+hap1]
    = local_marginals3[threadindex];
  }
  return;
}

__kernel void collect_xover_marginals_block_short(
__const int compressed_states,
__const int total_template_haplotypes,
__const int total_template_haplotype_pairs,
__constant int * block_marker,
__constant int * hap_pair,
__global packedshort_t * short_block,
__global float * backward_xover_marginals_big,
__local packedshort_t * local_short_prob_vec,
__local float * local_backward_vec,
__local float * local_marginals0,
__local float * local_marginals1,
__local float * local_marginals2,
__local float * local_marginals3
){
  int hap1 = get_group_id(0); 
  if (hap1>=total_template_haplotypes) return;
  int threadindex = get_local_id(0);
  local_marginals0[threadindex] = 0;
  local_marginals1[threadindex] = 0;
  local_marginals2[threadindex] = 0;
  local_marginals3[threadindex] = 0;
  barrier(CLK_LOCAL_MEM_FENCE);
  for(int chunk=0;chunk<=total_template_haplotypes/BLOCK_WIDTH;++chunk){
    int hapoffset = chunk*BLOCK_WIDTH;
    int hap2 = hapoffset+threadindex;
    if (hap2<total_template_haplotypes){
      int chunkwidth = total_template_haplotypes-hapoffset<BLOCK_WIDTH?total_template_haplotypes-hapoffset:BLOCK_WIDTH;
      if (threadindex<chunkwidth/2){
        local_short_prob_vec[threadindex] =  short_block[block_marker[0]*compressed_states+hap1*total_template_haplotype_pairs+hapoffset/2+threadindex];
      }
      barrier(CLK_LOCAL_MEM_FENCE);
      float newprob = short2float(local_short_prob_vec[threadindex/2].compressed[threadindex%2]);
      local_marginals0[threadindex] += (hap_pair[HAP_FIRST]==hap1 &&
      hap_pair[HAP_SECOND]==hap2)*newprob;
      local_marginals1[threadindex] += (hap_pair[HAP_FIRST]==hap1) *newprob;
      local_marginals2[threadindex] += (hap_pair[HAP_SECOND]==hap2) *newprob;
      local_marginals3[threadindex] += (hap1<=hap2)* newprob;
    }
    barrier(CLK_LOCAL_MEM_FENCE);
  }
  for(int s=BLOCK_WIDTH/2; s>0; s>>=1) {
    if (threadindex<s) {
      local_marginals0[threadindex]+=local_marginals0[threadindex+s];
      local_marginals1[threadindex]+=local_marginals1[threadindex+s];
      local_marginals2[threadindex]+=local_marginals2[threadindex+s];
      local_marginals3[threadindex]+=local_marginals3[threadindex+s];
    }
    barrier(CLK_LOCAL_MEM_FENCE);
  }
  if(threadindex==0){
    backward_xover_marginals_big[0*total_template_haplotypes+hap1]
    = local_marginals0[threadindex];
    backward_xover_marginals_big[1*total_template_haplotypes+hap1]
    = local_marginals1[threadindex];
    backward_xover_marginals_big[2*total_template_haplotypes+hap1]
    = local_marginals2[threadindex];
    backward_xover_marginals_big[3*total_template_haplotypes+hap1]
    = local_marginals3[threadindex];
  }
  return;
}

__kernel void collect_xover_marginals2_block(
__const int total_template_haplotypes,
__global float * xover_marginals_big,
__global float * xover_marginals,
__local float * local_marginals
){
  int marginal_type = get_group_id(1);
  int threadindex = get_local_id(0);
  local_marginals[threadindex] = 0;
  barrier(CLK_LOCAL_MEM_FENCE);
  for(int chunk=0;chunk<=total_template_haplotypes/BLOCK_WIDTH;++chunk){
    int hapindex = chunk*BLOCK_WIDTH+threadindex;
    if (hapindex<total_template_haplotypes){
      local_marginals[threadindex] += xover_marginals_big[
      marginal_type*total_template_haplotypes+hapindex];
    }
    barrier(CLK_LOCAL_MEM_FENCE);
  }
  for(int s=BLOCK_WIDTH/2; s>0; s>>=1) {
    if (threadindex<s) {
      local_marginals[threadindex]+=local_marginals[threadindex+s];
    }
    barrier(CLK_LOCAL_MEM_FENCE);
  }
  if( threadindex==0){
    xover_marginals[marginal_type] = local_marginals[threadindex];
  }
  return;
}

